#include <iostream>
#include <string>
#include <math.h>
#include <fstream>
#include "TFloat.h"
#include <stdio.h>
#include <stdlib.h>
#include <iomanip>

using namespace std;

int cantBis = 0;
int cantNew = 0;
int cantBisNew = 0;

int DEBUG = 1;
double GRAVEDAD = 9.81;

#define CantDeRebotes 2

/************************************************************************************************************/
/*                                               Sin Rozamiento                                             */
/************************************************************************************************************/

TFloat posicionSinRozamiento(TFloat h, TFloat vi, TFloat x, TFloat g, size_t precision) //Funcion posicion. Caso sin rozamiento
{
	TFloat DOS (2, precision);
    return (h + vi*x - (g/DOS)*x*x);
}

TFloat velocidadSinRozamiento(TFloat vi, TFloat x, TFloat g) //Funcion velocidad. Caso sin rozamiento
{
    return (vi - g*x);
}

TFloat veloDerivSin(TFloat vi, TFloat A, TFloat x, TFloat g, size_t precision) //Derivada de la velocidad. Caso sin rozamiento
{
	TFloat MENOS_UNO (-1, precision);
    return (MENOS_UNO*g);
}

TFloat impactoSin(TFloat fr, TFloat vi, TFloat x, TFloat g, size_t precision)
{
	TFloat MENOS_UNO (-1, precision);
    return (MENOS_UNO * fr * velocidadSinRozamiento(vi, x, g));
}

TFloat posSinRozamientoBiseccion(TFloat h, TFloat vi, TFloat a, TFloat b, TFloat ErrorBiseccion, TFloat g, size_t precision, int n)
{
    if ((posicionSinRozamiento(h, vi, a, g, precision) * posicionSinRozamiento(h, vi, b, g, precision)).dbl() >= 0)
    {
        cout << "FAIL (Posicion Sin Ronzamiento Biseccion)" << endl;
        cout << "f(a): " << posicionSinRozamiento(h, vi, a, g, precision).dbl() << endl;
        cout << "f(b): " << posicionSinRozamiento(h, vi, b, g, precision).dbl() << endl << endl;
        return 1;
    }

    TFloat m;
    int cant = 0;
    TFloat DOS(2, precision);

    while((b - a).dbl() > ErrorBiseccion.dbl() && cant < n)
    {
        m = (a + b)/DOS;
        TFloat x  = posicionSinRozamiento(h, vi, a, g, precision) * posicionSinRozamiento(h, vi, m, g, precision);
        if (x.dbl() < 0)
        {
            b = m;
        }
        else if (x.dbl() > 0)
        {
            a = m;
        }
        else
        {
            break;
        }

        cant++;
    }

    cantBis = cant;
    return m;
}

TFloat posSinRozamientoNewton(TFloat h, TFloat vi, TFloat xf, TFloat ErrorNewton, TFloat g, size_t precision)
{
    TFloat xi;
    int cant = 0;

    do
    {
        xi = xf;
        xf = xi - posicionSinRozamiento(h, vi, xi, g, precision)/velocidadSinRozamiento(vi, xi, g);
        cant++;
    }
    while(fabs((xf - xi).dbl()) > ErrorNewton.dbl());
    cantNew = cant;

    return xf;
}

TFloat velocidadBiseccionSin(TFloat vi, TFloat a, TFloat b, TFloat ErrorBiseccion, TFloat g, size_t precision, int n)
{
    if ((velocidadSinRozamiento(vi, a, g) * velocidadSinRozamiento(vi, b, g)).dbl() >= 0)
    {
        cout << "FAIL (Velocidad Biseccion Sin Rozamiento)" << endl;
        cout << "f(a) = "<< velocidadSinRozamiento(vi, a, g).dbl()<< endl;
        cout << "f(b) = "<< velocidadSinRozamiento(vi, b, g).dbl()<< endl;
        return 1;
    }

    TFloat m;
    int cant = 0;
    TFloat DOS(2, precision);

    while((b - a).dbl() > ErrorBiseccion.dbl() && cant<n)
    {
        m = (a + b)/DOS;
        TFloat x = velocidadSinRozamiento(vi, a, g) * velocidadSinRozamiento(vi, m, g);
        if (x.dbl() < 0)
        {
            b = m;
        }
        else if (x.dbl() > 0)
        {
            a = m;
        }
        else
        {
            break;
        }

        cant++;
    }

    cantBis = cant;
    return m;
}

TFloat velocidadNewtonSin(TFloat vi, TFloat A, TFloat xf, TFloat ErrorNewton, TFloat g, size_t precision)
{
    TFloat xi;
    int cant = 0;
    do
    {
        xi = xf;
        xf = xi - (velocidadSinRozamiento(vi, xi, g)/veloDerivSin(vi, A, xi, g, precision));
        cant++;
    }
    while(fabs((xf - xi).dbl()) > ErrorNewton.dbl());
    cantNew = cant;

    return xf;
}

TFloat eMecanicaSin(TFloat h, TFloat vi, TFloat t, TFloat g, size_t precision)
{
    TFloat DOS (2, precision);
    TFloat y = velocidadSinRozamiento(vi, t, g);
    return(g*posicionSinRozamiento(h, vi, t, g, precision) + (y * y/DOS));
}

/************************************************************************************************************/
/*                                               Con Rozamiento                                             */
/************************************************************************************************************/

TFloat posicionConRozamiento(TFloat h, TFloat vi, TFloat A, TFloat x, TFloat g, size_t precision) //Funcion posicion. Caso con rozamiento
{
	TFloat MENOS_UNO (-1, precision);
    return h + vi/A + g/(A*A) - (g*x)/A - (vi/A + g/(A*A))*(MENOS_UNO*A*x).exponencial();
}

TFloat velocidadConRozamiento(TFloat vi, TFloat A, TFloat x, TFloat g, size_t precision) //Funcion velocidad. Caso con rozamiento
{
	TFloat MENOS_UNO (-1, precision);
    return ((vi + g/A)*(MENOS_UNO*A*x).exponencial()) - g/A;
}

TFloat veloDerivCon(TFloat vi, TFloat A, TFloat x, TFloat g, size_t precision) //Derivada de la velocidad. Caso con rozamiento
{
	TFloat MENOS_UNO (-1, precision);
    return (vi+g/A)*(MENOS_UNO*A*x).exponencial()*MENOS_UNO*A;
}

TFloat impactoCon(TFloat fr, TFloat vi, TFloat A, TFloat x, TFloat g, size_t precision)
{
	TFloat MENOS_UNO (-1, precision);
    return (MENOS_UNO * fr * velocidadConRozamiento(vi, A, x, g, precision));
}

TFloat posConRozamientoBiseccion(TFloat h, TFloat vi, TFloat A, TFloat a, TFloat b, TFloat ErrorBiseccion, TFloat g, size_t precision, int n)
{
    if ((posicionConRozamiento(h, vi, A, a, g, precision) * posicionConRozamiento(h, vi, A, b, g, precision)).dbl() >= 0)
    {
        cout << "FAIL (Posicion Con Ronzamiento Biseccion)" << endl;
        cout << "f(a): " << posicionConRozamiento(h, vi, A, a, g, precision).dbl() << endl;
        cout << "f(b): " << posicionConRozamiento(h, vi, A, b, g, precision).dbl() << endl << endl;
        return 1;
    }

    TFloat m;
    int cant = 0;
    TFloat DOS(2, precision);

    while((b - a).dbl() > ErrorBiseccion.dbl() && cant < n)
    {
        m = (a + b)/ DOS;
        TFloat x  = posicionConRozamiento(h, vi, A, a, g, precision) * posicionConRozamiento(h, vi, A, m, g, precision);

        if (x.dbl() < 0)
        {
            b = m;
        }
        else if (x.dbl() > 0)
        {
            a = m;
        }
        else
        {
            break;
        }

        cant++;
    }

	cantBis = cant;
    return m;
}

TFloat posConRozamientoNewton(TFloat h, TFloat vi, TFloat A, TFloat xf, TFloat ErrorNewton, TFloat g, size_t precision, int n)
{
    TFloat xi;
    int cant = 0;
    do
    {
        xi = xf;
        xf = xi - posicionConRozamiento(h, vi, A, xi, g, precision)/velocidadConRozamiento(vi, A, xi, g, precision);
        cant++;
    }
    while(fabs((xf - xi).dbl()) > ErrorNewton.dbl()  && cant < n);
    cantNew = cant;

    return xf;
}

TFloat velocidadBiseccionCon(TFloat vi, TFloat A, TFloat a, TFloat b, TFloat ErrorBiseccion, TFloat g, size_t precision, int n)
{
    if ((velocidadConRozamiento(vi, A, a, g, precision) * velocidadConRozamiento(vi, A, b, g, precision)).dbl() >= 0)
    {
        cout << "FAIL (Velocidad Biseccion Con Rozamiento) " << endl;
        cout << "f(a): " << velocidadConRozamiento(vi, A, a, g, precision).dbl() << endl;
        cout << "f(b): " << velocidadConRozamiento(vi, A, b, g, precision).dbl() << endl;
        cout << endl;
        return 1;
    }

    TFloat m;
    int cant = 0;
    TFloat DOS(2, precision);

    while((b - a).dbl() > ErrorBiseccion.dbl() && cant<n)
    {
        m = (a + b)/DOS;
        TFloat x = velocidadConRozamiento(vi, A, a, g, precision) * velocidadConRozamiento(vi, A, m, g, precision);
        if (x.dbl() < 0)
        {
            b = m;
        }
        else if (x.dbl() > 0)
        {
            a = m;
        }
        else
        {
        	break;
        }

        cant++;
    }

    cantBis = cant;
    return m;
}

TFloat velocidadNewtonCon(TFloat vi, TFloat A, TFloat xf, TFloat ErrorNewton, TFloat g, size_t precision)
{
    TFloat xi;
    int cant = 0;

    do
    {
        xi = xf;
        xf = xi - (velocidadConRozamiento(vi, A, xi, g, precision)/veloDerivCon(vi, A, xi, g, precision));
        cant++;
    }
    while(fabs((xf - xi).dbl()) > ErrorNewton.dbl());
    cantNew = cant;

    return xf;
}

TFloat eMecanicaCon(TFloat h, TFloat vi, TFloat A, TFloat t, TFloat g, size_t precision)
{
    TFloat DOS (2, precision);
    TFloat y = velocidadConRozamiento(vi, A, t, g, precision);
    return(g*posicionConRozamiento(h, vi, A, t, g, precision)+(y * y/DOS));
}

/************************************************************************************************************/
/*                                               Tests                                                      */
/************************************************************************************************************/

void testBiseccionConRozamiento(TFloat h, TFloat vi, TFloat A, TFloat ErrorBiseccion, TFloat g, size_t precision, char *fileName)
{
    if (DEBUG) cout << "Entro a testBiseccionConRozamiento" << endl;

	ofstream myfile;
	myfile.open(fileName, ios::out | ios::app);

	int n = 10000;
	TFloat a (10, precision);
	TFloat b (150, precision);

	TFloat res = posConRozamientoBiseccion(h, vi, A, a, b, ErrorBiseccion, g, precision, n);
	myfile << "precision: " << precision << " t: " << setprecision(precision)  << res.dbl() << " alpha: " << A.dbl() << endl;

	myfile.close();

    if (DEBUG) cout << "Salgo de testBiseccionConRozamiento" << endl;
}

void testBiseccionSinRozamiento(TFloat h, TFloat vi, TFloat ErrorBiseccion, TFloat g, size_t precision, char *fileName)
{
    if (DEBUG) cout << "Entro a testBiseccionSinRozamiento" << endl;

	ofstream myfile;
	myfile.open(fileName, ios::out | ios::app);

	int n = 10000;
	TFloat a (0, precision);
	TFloat b (150, precision);

	TFloat res = posSinRozamientoBiseccion(h, vi, a, b, ErrorBiseccion, g, precision, n);
	myfile << "precision: " << precision << " " << "t: " << setprecision(precision)  << res.dbl() << " it: " << cantBis << endl;

	myfile.close();

    if (DEBUG) cout << "Salgo de testBiseccionSinRozamiento" << endl;
}

void testNewtonConRozamiento(TFloat h, TFloat vi, TFloat A, TFloat ErrorNewton, TFloat g, size_t precision, char *fileName)
{
    if (DEBUG) cout << "Entro a testNewtonConRozamiento" << endl;

	ofstream myfile;
	myfile.open(fileName, ios::out | ios::app);

	int n = 10000;
	TFloat x (150, precision);

	TFloat res = posConRozamientoNewton(h, vi, A, x, ErrorNewton, g, precision, n);
	myfile << "precision: "<< precision << " t: " << setprecision(precision)  << res.dbl() << " alpha: "<< A.dbl() << endl;

	myfile.close();

    if (DEBUG) cout << "Salgo de testNewtonConRozamiento" << endl;
}

void testNewtonSinRozamiento(TFloat h, TFloat vi, TFloat ErrorNewton, TFloat g, size_t precision, char *fileName)
{
    if (DEBUG) cout << "Entro a testNewtonSinRozamiento" << endl;

	ofstream myfile;
	myfile.open(fileName, ios::out | ios::app);

	TFloat x (150, precision);

	TFloat res = posSinRozamientoNewton(h, vi, x, ErrorNewton, g, precision);
	myfile << "precision: " << precision << " t: " << setprecision(precision)  << res.dbl() << " it: " << cantNew << endl;

	myfile.close();

    if (DEBUG) cout << "Salgo de testNewtonSinRozamiento" << endl;
}

void testBiseccionNewtonConRozamiento(TFloat h, TFloat vi, TFloat A, TFloat ErrorNewton, TFloat ErrorBiseccion, TFloat g, size_t precision, char *fileName)
{
    if (DEBUG) cout << "Entro a testBiseccionNewtonConRozamiento" << endl;

	ofstream myfile;
	myfile.open(fileName, ios::out | ios::app);

	int n = 10000;
	TFloat a (0, precision);
	TFloat b (150, precision);

	TFloat res = posConRozamientoBiseccion(h, vi, A, a, b, ErrorBiseccion, g, precision, n);
	res = posConRozamientoNewton(h, vi, A, res, ErrorNewton, g, precision, n);
    myfile << "precision: "<< precision << " t: " << setprecision(precision)  << res.dbl() << " alpha: "<< A.dbl() << endl;

	myfile.close();

    if (DEBUG) cout << "Salgo de testBiseccionNewtonConRozamiento" << endl;
}

void testBiseccionNewtonSinRozamiento(TFloat h, TFloat vi, TFloat ErrorNewton, TFloat ErrorBiseccion, TFloat g, size_t precision, char *fileName)
{
    if (DEBUG) cout << "Entro a testBiseccionNewtonSinRozamiento" << endl;

	ofstream myfile;
	myfile.open(fileName, ios::out | ios::app);

	int n = 10000;
	TFloat a (0, precision);
	TFloat b (150, precision);

	TFloat res = posSinRozamientoBiseccion(h, vi, a, b, ErrorBiseccion, g, precision, n);
    cantBisNew = cantBis;

	res = posSinRozamientoNewton(h, vi, res, ErrorNewton, g, precision);
    cantBisNew += cantNew;

	myfile << "precision: " << precision << " t: " << setprecision(precision)  << res.dbl() << " it: " << cantBisNew << endl;

	myfile.close();

    if (DEBUG) cout << "Salgo de testBiseccionNewtonSinRozamiento" << endl;
}

void testAlturaMaxima(TFloat h, TFloat vi, TFloat A, TFloat fr, TFloat ErrorNewton, TFloat ErrorBiseccion, TFloat g, size_t precision, char *fileName)
{
    if (DEBUG) cout << "Entro a testAlturaMaxima" << endl;

	ofstream myfile;
	myfile.open(fileName, ios::out | ios::app);

    TFloat primerRebote, segundoRebote, tiempoAlturaMaxima, alturaMaxima;
	int n = 10000;
	TFloat a (0.00001, precision);
	TFloat b (150, precision);

	// Primer rebote.
    primerRebote = posConRozamientoBiseccion(h, vi, A, a, b, ErrorBiseccion, g, precision, n);
    primerRebote = posConRozamientoNewton(h, vi, A, primerRebote, ErrorNewton, g, precision, n);
    primerRebote = primerRebote;

    // Segundo rebote y obtengo altura maxima.
    vi = impactoCon(fr, vi, A, primerRebote, g, precision); ///NO ENTIENDO POR QUE Vi DESPUES DEL IMPACTO ES NEGATIVO!!
    h = 0;
    segundoRebote = posConRozamientoNewton(h, vi, A, primerRebote, ErrorNewton, g, precision, 10000);

	segundoRebote = primerRebote + segundoRebote;
	tiempoAlturaMaxima = velocidadBiseccionCon(vi, A, a, b, ErrorBiseccion, g, precision, n);
	tiempoAlturaMaxima = velocidadNewtonCon(vi, A, tiempoAlturaMaxima, ErrorNewton, g, precision);
	alturaMaxima = posicionConRozamiento(h, vi, A, tiempoAlturaMaxima, g, precision);

	myfile << "precision: " << precision << " altura: " << alturaMaxima.dbl() << " 1erRebote: " << primerRebote.dbl() << " 2doRebote: " << segundoRebote.dbl() << " alpha: " << A.dbl() << " fr: " << fr.dbl() <<  endl;

    myfile.close();

    if (DEBUG) cout << "Salgo de testAlturaMaxima" << endl;
}

void testEnergiaMecanica(TFloat h, TFloat vi, TFloat A, TFloat fr, TFloat ErrorNewton, TFloat ErrorBiseccion, TFloat g, size_t precision, char *fileName)
{
    if (DEBUG) cout << "Entro a testEnergiaMecanica" << endl;

	ofstream myfile;
	myfile.open(fileName, ios::out | ios::app);

    TFloat primerRebote, segundoRebote, energia;
	int n = 10000;
	TFloat a (0.00001, precision);
	TFloat b (150, precision);

	// Primer rebote.
    primerRebote = posConRozamientoBiseccion(h, vi, A, a, b, ErrorBiseccion, g, precision, n);
    primerRebote = posConRozamientoNewton(h, vi, A, primerRebote, ErrorNewton, g, precision, n);
    primerRebote = primerRebote;
    energia = eMecanicaCon(h, vi, A, 0, g, precision);

    // Segundo rebote y obtengo energia.
    vi = impactoCon(fr, vi, A, primerRebote, g, precision); ///NO ENTIENDO POR QUE Vi DESPUES DEL IMPACTO ES NEGATIVO!!
    h = 0;
    segundoRebote = posConRozamientoNewton(h, vi, A, primerRebote, ErrorNewton, g, precision, 10000);
    energia = energia - eMecanicaCon(h, vi, A, primerRebote, g, precision);

	segundoRebote = primerRebote + segundoRebote;

	myfile << "alpha: " << A.dbl() << " fr: "<< fr.dbl() << " precision: "<< precision << " energia: " << energia.dbl() << endl;

    myfile.close();

    if (DEBUG) cout << "Salgo de testEnergiaMecanica" << endl;
}

/************************************************************************************************************/
/*                                               Main                                                       */
/************************************************************************************************************/

int main(int argc, char *argv[])
{
    char *fileName;

	if (argc == 11)
	{
		int test = atoi(argv[1]);
		size_t precision = (size_t)atof(argv[2]);
		TFloat h ((double)atof(argv[3]), precision);
		TFloat vi ((double)atof(argv[4]), precision);
		TFloat m ((double)atof(argv[5]), precision);
		TFloat cr ((double)atof(argv[6]), precision);
		TFloat A = cr/m;
		TFloat fr ((double)atof(argv[7]), precision);
		TFloat ErrorNewton ((double)atof(argv[8]), precision);
		TFloat ErrorBiseccion ((double)atof(argv[9]),precision);
		TFloat g (GRAVEDAD, precision);
		fileName = argv[10];

		cout << "test: " << test << " h: " << h.dbl() << " v0: " << vi.dbl() << " m: " << m.dbl() << " cr: " << cr.dbl();
		cout << " fr: " << fr.dbl() << " alpha: " << A.dbl() << " ErrorNewton: " << ErrorNewton.dbl();
		cout << " ErrorBiseccion: " << ErrorBiseccion.dbl() << " precision: " << precision << " file: " << fileName << endl;

		switch (test) {
		case 1: // Biseccion.
			if (A.dbl() > 0)
				testBiseccionConRozamiento(h, vi, A, ErrorBiseccion, g, precision, fileName);
			else
				testBiseccionSinRozamiento(h, vi, ErrorBiseccion, g, precision, fileName);
			break;
		case 2: // Newton.
			if (A.dbl() > 0)
				testNewtonConRozamiento(h, vi, A, ErrorNewton, g, precision, fileName);
			else
				testNewtonSinRozamiento(h, vi, ErrorNewton, g, precision, fileName);
			break;
		case 3: // Biseccion y Newton.
			if (A.dbl() > 0)
				testBiseccionNewtonConRozamiento(h, vi, A, ErrorNewton, ErrorBiseccion, g, precision, fileName);
			else
				testBiseccionNewtonSinRozamiento(h, vi, ErrorNewton, ErrorBiseccion, g, precision, fileName);
			break;
		case 4: // Altura maxima luego del primer y segundo rebote.
			testAlturaMaxima(h, vi, A, fr, ErrorNewton, ErrorBiseccion, g, precision, fileName);
			break;
		case 5: // Energia Mecanica.
			testEnergiaMecanica(h, vi, A, fr, ErrorNewton, ErrorBiseccion, g, precision, fileName);
			break;
		default: // Si no matchea con ningun test.
			cout << "No hay ningun test asociado al numero pasado por parametro" << endl;
			break;
		}
	}
	else
	{
		cout << "Te faltan parametros" << endl;
	}

    return 0;
}
